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(57) Abstract: A method for estimating near-surface material properties in the 
vicinity of a locally dense group of seismic receivers is disclosed. The method 
includes receiving seismic data that has been measured by a locally dense group 
of seismic receivers. Local derivatives of the waveficld are estimated such that 
the derivatives are centered at a single location preferably using the l.ax-Wen- 
droff correction. Physical relationships between the estimated derivatives in- 
cluding the free surface condition and wave equations are used to estimate near- 
surface material properties in the vicinity of the receiver group. Another em- 
bodiment ol' the invention is disclosed wherein the physical relationships used 
to estimate material properties are derived from the physics of plane waves arriv- 
ing at the receiver group, and the group of receivers does not include any buried 
receivers. 
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System and Method for Estimating Seismic Material 

Properties 

FIELD OF THE INVENTION : 
5 The present invention relates, to the field of 

determining characteristics of the land seismic 
wavefield. In particular, the invention relates to a 
system and method of processing seismic data from 
locally dense groups of receivers to estimate seismic 
10 material properties close to the receivers. 

BACKGROUND OF THE INVENTION : 

Previous studies have developed and tested 
methods to invert aspects of seismic data for elastic 

15 properties of Earth structure that are located either 
close to, or far from the seismic data receivers- When 
the data inverted are the actual seismic waveforms 
recorded, this process is called waveform inversion. 
These techniques compare seismic waveforms generated 

20 for a . synthetic Earth model with those recorded; the 
model is then updated in such a way that the misfit 
between the two sets of waveforms is reduced. Hence, up 
to now waveform inversion studies have solved the wave 
equation in an Earth model to predict and then fit the 

25 particular waveforms that were recorded, where the 

predicted waveforms are derived from the wave equation 
solution. Such techniques include ground roll inversion 
and diving wave tomography. However, these techniques 
suffer from at least one of the following shortcomings: 

30 non-uniqueness of solution, sensitivity to noise, and 
unpractical source receiver configurations. 



WO 01/53853 



- 2 - 



PCT/GB01/00182 



SUMMARY OF THE INVENTION ; 

• Thus, it is an object of the present 
invention to provide a method of processing seismic 
data from locally dense groups of receivers to estimate 
5 seismic material properties close to the receivers. 
According to the invention, a method of 
estimating near-surface material properties in the 
vicinity of a locally dense group of seismic receivers 
is provided. The method includes receiving seismic 
10 data that has been measured by a locally dense group of 
seismic receivers, and the data " representing earth 
motion caused by the seismic wavefield. Local 
derivatives of the wavefield are estimated such that 
the derivatives are centered at a single location in 
15 the vicinity of the receiver group. Physical 

relationships between the estimated derivatives are 
used to estimate near-surface material properties in 
the vicinity of the receiver group. 

According to a preferred embodiment of the 
20 invention, the derivatives are centered using a 

technique substantially similar to the Lax-Wendrof f 
correction, and the physical relationships used to 
estimate material properties include the free surface 
condition. Also according to a preferred embodiment, 
25 the physical relationships used to estimate material 

properties include wave equations, and at least one of 
the receivers in the receiver group is buried below the 
earth surface such that the receiver group encloses a 
volume . 

30 According to an embodiment of the invention, 

the receivers in the locally dense group are spaced 
apart by around 1 meter or less, or according to 
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another embodiment/ by about one-fifth of the shortest 
wavelength of interest or less. 

According to another embodiment of the 
invention, the receivers are located at or near the sea 
5 bottom, and the physical relationships used to estimate 
material properties include fluid-solid boundary 
conditions. 

According to another preferred embodiment of 
the invention, the physical relationships used to 
10 estimate material properties are derived from the 

physics of plane waves arriving at the receiver group, 
and the group of receivers does not include any buried 
receivers . 

The invention is also embodied in an 
15 apparatus for estimating near-surface- material 

properties in the vicinity of a locally dense group of 
seismic receivers. 

BRIEF DESCRIPTION OF THE DRAWINGS : 
20 Figure 1 shows an isotropic Earth model used 

in the description of the invention; . 

Figures 2a-d illustrate various geometries 
for the receiver groups, according to the invention ,- 

Figure 3 is a schematic illustration of a 
25 seismic data acquisition and processing system, 

according a preferred embodiment of the invention; 

Figures 4a-b illustrates the effect that mis- 
centering has on vertical derivative estimates; 

Figure 5 illustrates the frequency dependency 
30 of the contribution of the Lax-Wendrof f correction, 

according to a preferred embodiment of the invention; 

Figure 6 shows the misfit function jE7( 35 ) for 
the medium and complete time series data of which 
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snapshots are illustrated in figures 4a and 4b, 
according to the invention; 

Figure 7 is a flow chart showing the steps of 
estimating material properties according to preferred 
5 embodiments of the invention; and 

Figure 8 is a flow chart showing further . 
detail of the preferred methods for estimating material 
properties, according to preferred embodiments of the 
invention. 

10 

DETAILED DESCRIPTION OF THE INVENTION : 

Volumetric recording of the wavefield can be 
used to calculate and to invert spatial wavefield 

15 derivatives for P and S velocities in the Earth in the 
neighbourhood of a small, closely-spaced array of 
receivers. "Volumetric recording" refers to the fact 
that the array approximately encloses a volume of the 
Earth. This allows spatial derivatives of the wavefield 

20 in all directions to be calculated. The Quantities 

estimate<d are the effective velocities of the P and S 
components of the wavefield at any point in time. 
Hence, these will vary with both wave type and 
wavelength. If estimated for the near surface Earth 

25 structure, such velocities may be useful for statics 

estimation, or for separation of the wavefield into up- 
and down-going components. For further detail on 
separation of the wavefield see UK Patent Application 
GB 2 333 364 A (published 21 July 1999), and UK Patent 

30 Application entitled "System and Method for Seismic 
Wavefield Separation" (UK Patent Application No. 
0003406.6), filed on 15 February 2000, both of which 
are. incorporated herein by reference. 
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According to the invention, several inversion 
methods are provided. According to one method, 
referred to herein as "full wave equation inversion", 
all wavef ield derivatives that feature in the wave 
5 equation are measured explicitly so that the full wave 
equation can be inverted. This method can be applied 
anywhere within the Earth, but derivative measurement 
may be difficult logistically if the properties are 
required at depth beneath the surface of the Earth. 

10 According to the invention, the logistics may be 

simplified greatly if velocities are required in the 
near surface. . This simplification involves a 
transformation of the wave equation and introduces 
corrections, preferably analogous to the Lax-Wendroff 

15 corrections used in synthetic finite-difference 

techniques, to the wavef ield processing and inversion 
methodology. 

According to another embodiment of the 
invention, another inversion method is provided that is 

2 0 particularly applicable when estimating velocities 

close to the land surface. This method involves using 
the free surface conditions on the wavef ield to 
calculated first order vertical derivatives of the 
wavef ield that are consistent with the first order 
25 horizontal derivatives. These are then* equated with 
measured vertical derivatives giving a set of linear 
equations for P and S velocity (one equation from each 
point in time) These can be solved numerically as a 
function of time and frequency. 

3 0 To characterise sub~land reservoirs it is 

often necessary to analyse seismic data recorded on the 
land surface. Such data is contaminated with various 
forms of noise that must be removed before a correct 
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analysis can be carried out. Noise types specific to 
land-recording are derived from two main factors: 
first, static variations in the data -from receiver 
group to receiver group can be very large due to 
5 changes in sub-receiver group velocities and to 
variations in base topography of the uppermost 
weathered layer of the earth. Second, the land surface 
is exactly the point at which the up-going wavefield 
(including the signal from the reservoir) is reflecting 
10 and converting into a down-going wavefield. Data 

recorded at the earth's surface contains both up- and 
down-going wavefields; the up-going wavefield. must 
therefore be isolated in order to analyse the nature 
and true amplitudes of the signal coming up from the 
15 reservoir. This isolation is possible using techniques 
described in UK Patent Application entitled "System and 
Method for Seismic Wavefield Separation" (UK Patent 
Application No. 0003406.6), filed on 15 February 2000, 
but ordinarily requires sub-receiver group P and S 
20 velocities to be known. Hence, both types of noise can 
be attenuated if we are able to estimate sub- receiver 
group P and S velocities. 

Previous studies have developed and tested 
methods to invert aspects of seismic data for elastic 
25 properties of Earth structure that are located either 
close to, or far from the seismic data receivers.. When 
the data inverted are the actual seismic waveforms 
recorded, this process is referred to as waveform 
Inversion, These techniques compare seismic waveforms 
3 0 generated for a synthetic Earth model with those 

recorded; the model is then updated in such a way that 
the misfit between the two sets of waveforms is 
reduced. Hence, conventional waveform inversion 
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techniques have solved the wave equation in an Earth 
model to predict and then fit the particular waveforms 
that were recorded, where the predicted waveforms are 
derived from the wave equation solution. 
5 By contrast, according to the present 

invention, a method is provided that allows one to 
invert the complete wave equation for Earth properties 
rather than inverting any derived waveform quantity. 
• The elastic, isotropic wave equation for" particle 
10 displacement u (actually three dependent equations) can 
be written as: 

pu = f + (A + 2/z)V(V ■ u) - jjlV X (V x u), (1) 

15 where p is the density, X and ju are Lamp's 

constants, f denotes the distribution of body forces 
and V = [3i,32,d3] T where 3* is the derivative operator 
with respect to coordinate direction x*. Full wave 
equation inversion involves measuring all derivatives 

2 0 of the wavefield included in equation (1) ; once these 

and the body forces are known, equation (1) comprises a 
set of linear equations that may be solved for P and S 
. velocity a and P respectively where, 

When full wave equation inversion is 
performed at the Earth's surface (e.g., for land 
seismics) , estimates of third order horizontal 

3 0 derivatives of the wavefield are ordinarily required. 

According to a preferred embodiment of the invention, 



WO 01/53853 



- 8 - 



PCT/GBO 1/00 182 



another inversion method provides for the construction 
of a set of linear equations in a and |3 specifically 
for use in land seismics using only second order 
derivative estimates . 
5 The particular values of a and 0 that are 

recovered using either' method are those which control 
the propagation of waves past the receivers. Hence, 
these are effective P and S velocities averaged over a 
wavelength or so from the receivers, and may depend on 
10 both frequency and wave-type. 

Wave equation inversion is particularly 
suited to new acquisition geometries involving receiver 
arrays that allow spatial and temporal derivatives of 
the wavefield to be recorded explicitly. For example, 
15 refer to UK Patent Application No. 9921816.6,. filed 15 
September 1999, incorporated herein by reference* From 
studying the stability of such recording methods under 
perturbations in receiver locations, it is believed 
that whereas estimates of divergence and curl of the 
20 wavefield from volumetric recordings are relatively 
robust" under perturbations in element location or 
amplitude, they are more sensitive to perturbations in 
the orientation of the individual recording sensors.. 
Techniques based on these (and simpler) recording 
25 geometries can be used for separation of P and S waves 
and of upgoing and downgoing wavefields when recordings 
are, made on, or close to, the ground surface. If 
temporally distinct P wave arrivals can be 
distinguished on signals recorded in such a 
30 configuration at the Earth's surface then estimates of 
material properties close to the receivers can be 



WO 01/53853 



- 9 ~ 



PCT/GBO1/00182 



obtained using (assumed) plane, wave approximations of 
the incoming field - 

According to the present invention, these 
acquisition techniques can be adapted so that all 
5 wavefield derivatives pertinent to the wave equation, 
and hence to wave propagation local to the array, can 
be estimated. Using either preferred embodiment, these 
derivatives can be inverted to constrain effective 
material properties close to the receivers using 

10 recorded energy from any wave types, whether arrivals 
are temporally distinct or contemporaneous . These 
effective material properties are exactly those that 
can be used for separation of up- and down-going 
components of the wavefield. 

15 At various stages herein exemplary results 

using synthetic data will be presented. The data was 
generated using a plane-layered reflectivity code. See 
e.g., Kennett, B. L. , 1983, Seismic wave propagation in 
stratified media: Cambridge University Press,' 

2 0 Cambridge. The reflectivity code used works ' to an 

internal precision of 4 significant figures. Figure 1 
shows an isotropic Earth model used in the description 
of the invention. A Ricker wavelet with 50 Hz central 
frequency was injected at the source location. 
25 Figure 2 illustrates various geometries for 

the receiver groups, according to the invention. As 
the name suggests, volumetric wavefield recording 
involves recording aspects of the wavefield at points 
that approximately enclose a volume in the Earth. 

3 0 Figure 2(a) shows a geometry that approximately 

encloses a four point tetrahedral volume defined using 
the receiver locations as vertices. It has been shown 
that this geometry was sufficient to estimate the 
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average divergence, and curl of the wavefield within the 
tetrahedron, thus performing approximate separation of 
P and S energy (respectively divergence and curl of the 
wavefield) . See, Johan 0. A, Robertsson and Everhard 
5 Muyzert . Wavefield separation using a volume 

distribution of three component recordings. Geophys . 
Res. Lett., 26(18) :2821-2824 / 1999, incorporated herein 
by reference. 

Using the geometry shown in Figure 2(a), all 

10 first order spatial derivatives of the wavefield can be 
calculated simply by taking finite differences 
(hereinafter *fd's") between different components. For 
example, if the particle velocity v (i) is recorded at 
receiver i, and receivers 1, 2 and 3 in Figure 2(a) are 

15 located on the corners of an equilateral triangle with 
centroid vertically above receiver 4, then 



a 3 v „, 1 _ IU. ,3, 



0 where Cartesian coordinate directions are x lt 

x 2 and x 3 with x 3 vertically downwards, where in our 
notation d 3 v represents dv/dx^, and where Ax 3 is the 
depth of receiver 4. Other derivatives may be 
approximated using different fd stencils. 

5 Considering the aim to invert simultaneously 

various different wavefield derivatives for material 
properties using the wave ec~~atiori, one disadvantage of 
this geometry is immediately apparent. The depth 
derivative estimate in equation (3) is centered 

0 vertically above receiver 4 at depth dx^/2 whereas 
estimates of both 3iV and 9 2 v will use fd stencils 
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involving only receivers 1-3 and hence will be centered 
on the surface. The wave equation (1) relates 
derivatives at identical locations, and although this 
deviation in centering may seem small, it can be very 
5 significant as illustrated below. Ta obtain derivative 
estimates that are all centered at exactly the same 
location, at least two options are available: first, 
geophones might be planted in a more symmetric 
geometry. For example, equi-centered fd stencils for 

10 first-order spatial derivatives can easily be 

Constructed using the geometry shown in Figure 2 (b) . 
However, if such geometries are too labour intensive or 
are simply impractical to plant, a method that corrects 
the first-order derivative centering ' using higher-order 

15 derivatives as discussed in further detail below. 

To invert equation (1) for P and S velocity 
it is necessary to estimate second order spatial and 
temporal wavefield derivatives centered at identical 
locations. Computationally, the simplest way that this 

20 can be. achieved is to plant a receiver configuration 

similar to that in figure 2 (b) but this time using a 3 
x 3 array of receivers. Fd stencils similar as those 
discussed in further detail below can then be used to 
estimate second derivatives . 

25 In certain situations geometries using fewer 

receivers can be used (See, e.g., Figures 2(c) and 
2 (d) ) , but similarly to the first-order case above, if 
identically centered derivatives are to be measured 
then centering corrections should be made. -For the 

3 0 most part in this description, the simplified 

geometries are used in land seismic surveys where the 
free surface boundary condition allows certain 
derivatives to be calculated implicitly. However, 
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according to the invention, similar geometries can be 
applied in the ocean bottom or borehole seismics cases, 
where the fluid-solid boundary plus a pressure sensor 
allow the same derivatives to be calculated implicitly. 
5 A principal feature of the invention is that 

the receivers within each group used to calculate 
wavefield derivatives are spaced apart at most by" 
around one-fifth to one- tenth of the relevant 
wavelenth. Advantages of moving to smaller inter- 

10 receiver spacings include increased accuracy in the 
estimated derivatives, and greater validity in the 
assumption that the material properties are the same in 
the vicinity of the receivers. However, an advantage 
of wider spacing is less sensitivity to noise. These 

15 competing concerns in combination with the wavelength 
of interest (or more accurately, the projection of the 
wavefield on the recording surface) , should be taken 
into consideration when designing the receiver group 
spacing. According to a presently preferred 

2 0 embodiment, the locally dense receivers are spaced 

about 1 meter apart. However spacing can in some 
situations be larger, for e.g. around 2 meters, or 
smaller, e.g. 0.5 meters. According to a preferred 
embodiment the inter-receiver spacing is around 0.25 
25 meters or less. As mentioned, the optimal spacing of 
the receivers depends upon the wavelength of interest. 
There, should be at least two receivers within the 
projection of the shortest wavelength of interest on 
the recording surface. According to a preferred 

3 0 embodiment the receivers are spaced apart by a distance 

approximately equal to or less than one-fifth of the 
shortest wavelength of interest. 
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Figure 3 is a schematic illustration of a 
seismic data acquisition and processing system, 
according a preferred embodiment of the invention . 
Seismic sources 150, 152, and 154 are depicted which. 
5 impart vibrations into the earth at its surface 166. 
The vibrations imparted onto the surface 166 travel 
through the earth; this is schematically depicted in 
Figure 3 as arrows 17 0. The vibrations reflect off of 
certain subterranean surfaces, here depicted as surface 

10 162 and surface 164 , and eventually reach, and are 

detected by receiver groups 156, 158, and 160. Each 
receiver group comprises a number of receivers such as 
in the arrangements depicted in Figures 2a~d. 

Importantly, according to the preferred 

15 embodiment, the spacing of between the receivers within 
a single receiver group is substantially less than the 
spacing between the receiver groups. Schematically, 
this is shown in Figure 3 by the size 122 of receiver 
group 160 is substantially smaller than the distance 

20 120 between group 160' and an adjacent group 158. 

Each of the receivers in groups 15 6, 158, and 
160 convert the vibrations into electrical signals and 
transmit these signals to a central recording unit 170, 
usually located at* the local field site. Preferably, 

2 5 the data is not group formed, but data from each 

geophone component are recorded. The central recording 
unit typically has data processing capability such that 
it can perform a cross-correlation with the source 
signal thereby producing a signal having the recorded 

3 0 vibrations compressed into relatively narrow wavelets 

or pulses. In addition, central recording units may 
provide other processing which may be desirable for a 
particular application. Once the central processing 
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unit 170 performs the correlation and other desired 
processing, it typically stores the data in the form of 
time-domain traces on a magnetic tape. The data, in 
the form of magnetic tape is later sent for processing 
5 and analysis to a seismic data processing center, 

typically located in some other geographical location. 
The data transfer from the central recording unit 170 
in Figure 3 is depicted as arrow 17 6 to a data 
processor 180. 

10 Referring again to Figures 2a-d, receiver 

geometries similar to those in Figures 2(a) and 2(b) 
have the advantage that they may be planted at any 
orientation, and at any depth in the Earth, and will 
still provide the same derivative information . However, 

15 a disadvantage when deployed below the surface is that 
several receivers must be buried within the Earth 
(especially limiting if a 3- x 3 x 3 cubic geometry is 
used for second-order derivative estimation) . 
According to the invention, simpler receiver geometries 

2 0 can be used to obtain the same information when 
deployed at the Earth's surface. 

The elastic constitutive relation relates 
components of the stress tensor aft in a source free 
region to the strain tensor components 

25 

<*ij = CijktEki , (4) 

where aju are the elastic stiffnesses. Index 
values 1 and 2 correspond to horizontal coordinates x 2 
30 and X2 whereas index value 3 corresponds to the 
vertical direction downwards x 3 ; Using the Voigt 
notation, equation (1) can be written as: 
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where C is the symmetric stiffness matrix 
with 21 independent components. If we assume that 
material properties in the near-surface environment are 
isotropic, the stiffness matrix takes the following 
form: 
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(6) 



where X and \i are the Lame constants. Strain 
Eij is related to particle velocity v* through 



15 



(7) 



where d t denotes a time derivative. The free- 
surface conditions on the wavefield that must hold at 
the Earth's surface can be written: 
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= 0 3 i = 1,... ,3. • (g) 

■Equation (8) states that all components of 
the stress tensor acting across the (horizontal) free 
surface of the Earth vanish. Substituting 
equations (5) , (6) and (7)' into (8) provides three 
constraints on particle velocities at the free surface: 

d z v 2 == - d 2 v z , (10J 

ds Vl = -0^. (11) 



Thus the vertical wavefield derivatives can 
be expressed in terms of horizontal derivatives. Note 
that the material properties only occur in equation 

15 (9), not in equations (10) and (11). 

In land acquisition, spatial distributions of 
3C receivers at the surface allow for the computation 
of horizontal spatial derivatives of particle 
velocities (or time-derivatives thereof if particle 

20 acceleration is recorded, etc.). Advantageously, 

according to the invention, invoking the free surface 
conditions (equations (9) -(H)) enables the 
calculatation of first order vertical derivatives using 
the same configuration of 3C receivers at the surface, 

25 provided that the ratio between the P and S velocities 
is known. 

In fact, if either of the receiver geometries 
in Figures 2(c) or (d) is used, the free surface 
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15 



condition allows all second order derivatives of the 
wavefield to be estimated as shown in further detail 
below. Note that applying fd stencils to the geometry 
in Figure 2 (d) also allows all third order purely 
horizontal derivatives to be estimated. 

According to the invention, the free surface 
conditions (equations (9) -{ID) can be used to rewrite 
the. wave equation (1) into a different form that is 
valid at the free surface in an isotropic medium. 
Provided all of the derivative estimates are centered 
at the same point on the free surface, this equation 
can be inverted for near-surface material properties. 
The wave equation comprises the constitutive 
relation (4) and the equation of motion: 



diV = -V-<7 
P 



P 



-+■ 0>>.oii -f c?auai 



H2) 
(13) 



Let us define the following terms : 



20 
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(14) 
(15) 



The isotropic constitutive relation (4) , the 
free surface conditions (9) -(11) and the equation of 
motion (13) yields the following set of free surface 
25 wave equations when no body forces are present: 
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d if Vi 



2^ 2 ) (Viva) 



(16] 

(17) 
(18) 



• where the vertical derivatives are taken 
downwards. As used herein, repeated subscripts on the 
derivative operator d denote multiple derivatives 
(e.g., 3 33 = a 3 a 3 ) . The final term on the right hand side 
of each of equations (16) -(18) contains only horizontal 
derivatives. The only depth derivatives in these 
equations are contained in the first term on the right 
of each equation, and these terms only involve the pure 
second order depth derivatives 8 33 v\ Hence, these 
equations form three constraints that can be inverted 
directly for material parameters a and 0, and no mixed 
derivatives in depth are required. The latter point is 
important since it suggests that we can measure all 
relevant derivatives using only a single buried 
geophone (e.g., using the geometries in Figures 2(c) 
and 2(d)) using the techniques described in further 
detail below. 

If the receiver geometry used does not allow 
derivatives to be centered at identical locations, 
material properties obtained by inserting derivative 
estimates into either form of the wave equation (1) or 
equations (16) to (13) will be inaccurate. In surface 
seismics this principally affects estimates of 
derivatives with respect to depth; logistically it is 
usually easy to center all horizontal derivatives at a 
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single location - on the surface, but the same centering 
is impossible to achieve directly for vertical 
derivative measurements . 

Figures 4a and 4b illustrate the effect that 
5 mis-centering has on vertical derivative estimates. 
From Figures 4a and 4b, it can be seen that for 
conventional seismic frequencies and typical near 
surface velocities, differences in centering of only 
12 . 5 cm can cause large deviations in the estimates 

10 obtained. The first order vertical derivative estimated 
using the free surface condition is calculated directly 
from the horizontal derivative estimates using 
equation (9), and hence is centered exactly on the free 
surface* The vertical derivative estimated explicitly 

15 using f-d stencils centered only 12.5 cm below the 

surface is shown to provide very poor approximations to 
the surface value in the example given (this is the 
case in many realistic examples tested). Hence, often 
it will be necessary to correct the derivative 

2 0 centering before wave equation inversion can be 
applied. 

Figures 4a and 4b estimate of 3 3 v 3 for a 
synthetic P wave arrival (Figure 4a) and surface wave 
arrival (Figure 4b) using the experiment geometry in 

2 5 Figure 1 with a receiver array similar to that in 

Figure 2(d), but with receivers spaced 25 cm apart. In 
each case, the solid curve is the derivative calculated 
using the free surface condition in equation (9) and 
hence is centered on the free surface. The thin dashed 

30 curve is d' 3 v 3 defined by equation (53), below, using a 
first-order fd stencil in depth, and hence is centered 
12.5 cm below the ground surface. The thick dashed 
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curve is derivative ^v 3 re-centered to the free surface 
using the Lax-Wendroff correction in equations (23) 
and (24), and almost entirely overlays the black curve. 
The horizontal axis is time in milliseconds (ms) v , the 
5 vertical axis is the derivative value in 1/ms. 

•At any point Xo on the Earth's surface a 
Taylor expansion of velocity V can be used to write 



10 



Ax 2 

v(x 0 + Ax 3 ) = v(x 0 ; + 5 3 v(x 0 ) + — i d^v(Ko) -f O(Asl) U*) 



.where Ax 3 = [0,0,Ax 3 ] T for any small depth 
increment Ax 3 . Hereinafter, the order of the derivative 
being estimated will be referred to as usual by first 
.5 order, second order, etc., and the accuracy of the 
Taylor approximation used for the estimates will be 
referred to as o(ax*) , o(Ax%) , etc. Rearranging this 
equation gives, 



0 L 



v(xp -h Ax 3 ) - y(xo) 



(20) 



Hence, if 3 33 v(x 0 ) can be estimated then the 
conventional 0(Ax 3 ) fd stencil for d 3 v(xo + Ax 3 /2) in 
equation (2 0) (the square-bracketed term, right hand 
side) can be modified to give 3 3 v(x 0 ) up to o(Ax 3 ). The 
correction term (second term, right hand side) consists 
of a linear (in Ax 3 ) approximation to the error 
introduced by mis-centering the derivative estimates 
using the conventional stencil. 
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Derivatives 3 33 v(x 0 ) can be related to 
horizontal and temporal derivatives using the wave 
equation. Since we wish to use these derivatives in a 
correction centered at the ground surface we will use 
5 the free surface wave equations (16) to (18) to obtain 
0 (Ax 3 ) corrections, Li say, defined to be equal to the 
second term on the right side of equation (20). These 
are obtained simply by rearranging equations (16) to 
(18) : 

2 ( 1 ~^) ai(V *' V * ) (21) 
2 (l-^) ^(V r v^ (22) 

«»«• = ^ + ( i - 2 ^) v » < * <23) 

Li = -^c^v (24) 

In modelling wave propagation using finite 
15 - difference techniques, according to the invention, the 
preferred method of computing this type of correction 
term using the wave equation to calculate any 
derivatives that can not be measured is known as Lax- 
Wendroff correction. For further detail on this method 
20 refer to: P. D. Lax and B. Wendroff, Difference schemes 
for hyperbolic equations with high order of accuracy. 
Comm. Pure appl . Math., 27, 1964, incorporated herein 
by reference See also, J. 0. Blanch and J. O. A. 
Roberts son, A modified Lax-Wendroff correction for wave 
25 propagation in media described by Zener elements. 



#33 1>2 
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Geophys. J\. Int., 131:381-386, 1997, also incorporated 
herein by reference. 

Using the Lax-Wendroff correction for 
correcting measured data advantageously provides 
5 increased ^accuracy in vertical wavef ield derivatives 
centered exactly on the free surface of the earth. 

Figures 4a and 4b show examples of the Lax- 
Wendrof f correction Li applied to first order 

derivative .estimates of 3aV3 using fd stencils centered 

10 12.5 cm below the surface. The corrected derivative 
estimate is virtually identical to that obtained 
exactly on the free surface (using the free surface 
conditions) . Hence, as long as the second order 
horizontal and temporal derivatives in equations (21) 

15 to (23) can be estimated with sufficient accuracy, all 
first order derivatives can be estimated to o(Ax*) 
exactly at the free surface. 

Figure 5 shows that the contribution of the 
Lax-Wendrof f correction for 8 3 v 3 is highly frequency - 

20 dependent/ In particular, it does not have a maximum 
amplitude at the same frequency as the maximum of the 
signal itself. This is because Li contains second 
derivatives in both space and time (the second order 
time derivative is equivalent to multiplication by -CO 2 

25. in the frequency domain) . 

Similarly to the case just described, the 
second order derivative estimates obtained using 
equation (54) (set out below) , are centered at a depth 
Ax 3 /4 below the surface . Estimates centered at the 

3 0 surface can be obtained by expanding the Taylor series 
in equation (20) to one additional term. Rearranging 
the terms yields: 



WO 01/53853 



- 23 - 



PCT/CBO 1/00182 



v(xo + A X< )-v(xo)_ 53v(xo) 



As 3 

gUv(x 0 ) ■+ O(Axl) (25) 



Hence, the 0(Ax 3 ) Lax-Wendrof f correction .Ii 2 
5 for second order derivatives is given by 



L 2 = — — — «5a33V (26} 



The first term on the right hand side of 
10 equation (25) is the approximation to d 33 v (x 0 + Ax 3 /4) 

given by equation (54) . Equation (26) provides a linear 
(in Ax 3 ) correction that shifts the derivative 
centering tox 0 - The third order derivatives with depth 
can be obtained by differentiating the wave 
15 equation (1) once with respect to x 3/ and again 

applying the free surface condition to simplify mixed 
derivatives in depth. Notice that we can not achieve 
this by differentiating expressions (21) to (23) with 
respect to depth since these expressions are only valid 
20 exactly on the free surface. The third order vertical 
derivatives at the free surface can be written: 
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^333^3 = 



^333^1 = 




(29) 



(27) 



(28) 



Hence we have expressed third order vertical 



derivatives in terms of second order time derivatives 
5 of first order horizontal derivatives {first terms on 
right of each equation) and third order horizontal 
derivatives (second term on right hand sides) . These 



isotropic P and S velocities gained using this 
preferred technique described herein is the maximum 
possible such information from any waveform inversion 
technique. This is because the derivative data (plus 

15 isotropic assumption) includes sufficient information 
to determine completely' the local wavefield within the 
recorded bandwidth. Hence/ there exists no other 
information about the wavefield that could possibly be 
used to improve constraints on the P and S velocities. 

2 0 It is possible to derive a set of linear 

constraints on the P and S velocities that involve only 
first order spatial wavefield derivatives. We do so by 
calculating the vertical derivatives of the wavefield 
both using, and without using the free surface 

25 condition in equation (9) . In each of these cases we 
obtain respectively: 



can be substituted into equation (26) to obtain the 
required Lax-Wendroff correction. 



10 



It is believed that the information on 
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(dat/a )fd 



s 




(33) 



(34) 



where fs and fd refer to "free surface" and 



5 



"finite difference" respectively, and v* is the x 3 



component of velocity at distance Ax 3 below the point 
at which v 3 is recorded. Equation (33) is really a 
definition of (<? 3 v 3 )fs using equation (9), and 
equation (34) is a finite difference approximation to 
10 the vertical derivative ^v 3 . Notice that in 

equation (34) , the third component of the correction Li 
has been applied to ensure that the derivative 
estimates are both centered at exactly the same 
location. 

15 If estimates of the derivatives on the right 

hand sides of equations (33) and (34) are made using fd 
stencils that are centered at approximately the same 
location, then two estimates of the same quantity are 
obtained, one of which uses an explicit depth 

20 derivative, one of which does not. Hence, 




(35) 



25 



provides a set of constraints on oe and p. 
Expanding equation (35) gives: 
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dr 3 

2 




(36) 



-J!? 2 [2(Vjy • y h ) - <k*(Vi©0] (37) 

where ^v 3 is given by equation (53) . This 
shows that a set of linear equations is obtained 
5 relating (3 2 and a? , one equation for each point in time. 
These can fc>e solved numerically. 

It will now be shown how equation (1) , 
equations (16) -(18) or equation (35) can be inverted to 
constrain P and S velocities, a and 3 respectively, 

10 given the derivative estimates. Equations ( 1 ) and ( 35 ) 

each result in a set of linear equations in a 2 and p 2 . 
These, systems can be solved very simply using standard 
techniques (e.g., W.. Menke, Geophysical data analysis: 
Discrete inverse theory (Revised edition) , volume 45 of 

15 International Geophysics Series. Academic Press Inc., 
Harcourt Brace Jovanovich, 1989; F. Press, B. P. 
Flannery, S. A. Teukolsky, and W. T. Vetterling. 
Numerical Recipies in FORTRAN, the art of scientific 
computing (2nd edition) . Cambridge University Press, 

20 1992; and R. L . Parker. Geophysical inverse theory. 

Princeton University Press, 1994.) However, in order 
to illustrate the solution uncertainty graphically we 
• will construct and display misfit functions in the time 
domain : 

25 
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Here I* N and R N refer respectively to the 
left- and right-hand sides of equation number (N) . 
5 These are vectors consisting of samples in time or 

frequency where the signals may have been filtered in 
time to extract any desired wavefield arrival, and 
filtered in frequency to extract any desired frequency 
component. The derivative estimates in each equation 

10 (N) are assumed to have been centered at identical 
locations. The factor (a 11 ) 2 represents the scalar 
variance of elements of L N , and its reciprocal is a 
suitable weighting factor so that contributions from 
signals with large amplitudes do not dominate E. 

15 Given sufficiently accurate derivative data, 

we can estimate parameters a and J3 by minimising any of 
J5? ( i), Sas), S ( i7), Ens), or £7(35)/ or a combination of all 
of these. Thus we ensure that any of equations (1) , 
(16) -(18) and (35) are satisfied as accurately as 

20 possible. 

First we construct a frequency domain misfit 
function based on the requirement that divergences 
satisfy equation (35) : 

f [|^)| -|j;(^IFpW|- |ii<«>|n ' 

25 s < J5 ' = l0S \ }■ (39) 

Here the modulus signs imply that only the 
Fourier domain amplitudes (within the signal frequency 
band) of left and right hand sides of equation (3 5) are 
3 0 compared. This removes problems associated with points 
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10 



in time where either v /IS or v yiQ in 

equation . (35 ) become numerically equivalent to zero. 

Figure 6 shows the misfit function E { i S) for 
the medium and complete time series data of which 
snapshots are illustrated in figures 4a and 4b. It is 
immediately apparent that the misfit function has a 
well defined minimum at the correct solution (a* = 
ISOOm/s, P b = 500m/s) . Hence, if sufficiently accurate 
data are available, the constraints offered by 
equation (35) (at each point in time) are sufficient to 
estimate both P and S velocities. 

The accuracy of such estimates will decrease 
with decreasing time series length or frequency 
content. Even using the complete modelled time series 
15 in this case there exists a linear trade-off that 
degrades the estimate of (3 more than that of a. 
However, noting that the colour scale is logarithmic, 
the. minimum in this surface still gives at least a 
factor of ten better fit to the data than oc t ±2% or 0 t 
2 0 ±7%. Hence, since time series from many sources can be 
inverted for the same P and S velocities, the estimate 
accuracy should be sufficient for most requirements. 

The preferred technique for estimating 1st 
4 and 2nd order spatial wavef ield derivatives will now be 
25 discussed in further detail. Using either of the 
receiver geometries in Figures 2(c) or 2(d) it is 
possible to estimate all first and second order 
derivatives in the wavef ield. Using receiver geometry d 
for instance, we may estimate horizontal derivatives 
30 centred at the free surface above the central point 
using the finite difference formulae: 
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c?i V 



1 

2 



24 Aan 1 
obtained by rotation of the above 



{ 

i 

2 



obtained from djV and using 

the free surface condition (e.g., equations (9)-(l 1)) 



Ta^ + 



obtained by rotation of the above 



.(11) _ v (io) 



,(7) _ V W 



Ar 2 



A^i 



+ 0(Aj:?,Ax 2 2 ) 



(40) 
(41) 
(42) 

(43) 

(44) 
(45) 



where bracketed superscripts denote the 
receiver number in figure 2(d) . Second mixed 
derivatives in the vertical direction can be obtained 
by using the free surface condition: 



(46) 
(47) 
(48) 
(49) 
(50) 



dl3Ul = 


-3llt>3 




d 23 t>;t = 


-<9l 2 t>3 




dl3"2 = 


-dl2t?3 




<?23 u 2 — 


— #22 1>3 




9iaVa = 






^23^3 = 


-d 2 (y H 





(51) 



10 All first and second order derivatives above 

can be estimated using only receivers on the surface. 
Finally, however, second order pure derivatives in 
depth must be obtained using both the free surface 
condition and receiver 17 at depth as follows. 
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5 



Define the velocity v s to be that estimated at the free 
surface directly above receiver 17, 

. V W _|_ v (7) + v (10) . y(ll) 

v J = ~ • (52) 

The finite-difference first order derivatives 
in depth (denoted d z ) are centred vertically above 
receiver. 17 half way to the surface: 

c%v = [v<"> - v'] J- + 0(Ax|) (S3) 



In addition, equation (42) estimates the same 
derivative at the surface. Hence, the difference 
between these two estimates can be used to obtain the 
15 second order depth derivatives centred at depth <£c 3 /4: 



20 



Thus, we may estimate all first and second 
order derivatives of the wavefield using the receiver 
geometry in Figure 2 (b) . For any other surface 
geometry (e.g., the minimal 5 -star spread in 
figure 2(c)), only the horizontal spatial derivative 
finite difference estimations (40) to (45) and the 
25 surface velocity estimate for V s in equation (52) 

change, whereas equations (46) to (51) and (53) to (54) 
remain the same. 

According to another embodiment of the 
invention, another method is provided in which spatial 
30 and temporal derivatives of the wavefield can be used 
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to estimate near-surface material properties. This 
method is also particularly applicable where dense 
recordings of the wavefield are available to allow 
spatial derivatives to be calculated and the wave types 
5 present in a specific temporal window to be known. 

Additionally, this method advantageously can be applied 
to data from surveys having the entire receiver group 
on the ground surface, i.e. where there is no buried 
receivers in the receiver group. 

10 According to this embodiment, the interaction 

of the wavefield with the free surface can easily be 
expressed in a potential formulation,. According to 
this embodiment, it is assumed that the selected 
temporal window contains one up-going plane P-wave and 

15 the positive z-direction is chosen to be downwards. The 
incoming P-wave generates both a reflected P- and a 
converted SV-wave at the surface, with corresponding 
scalar potentials : 

^ , [. /ski cost \1 

= A csp hw f - * - tj J (55) 

[ /sins , cos i \1 
I t + — — z - t 1 (56) 

where a and P are the P- and S-velocity respectively, i 
is the angle of incidence of the incoming P-wave and j 
the angle of reflection for the converted S-wave. The 
25 corresponding expressions for the particle velocities 
can be obtained from equations (55-56) through: 
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_ $&inc $&r*fl _ d^refl ( 58 ) 

= £|i!!l + +^re fi (59) 

nix fix j9rr 



If ratios between u* and u z or derivatives of 
these are calculated the exponential expressions will 
vanish for recordings at the free surface, where z ~ 0 
This yields : 



tt* (Vp.p — 1) cost + sinV V>5 

d*a s fdx cos t sim (1 — V>j?) — Vps 

d**}dt ~ a ^ f ^ + 7w) __ ./^-^ (61) 



where y = a/|3 and V PP and Vps are the reflection 
and conversion coefficients for the potentials at the 
free surface. Equations (60), (61) and (62) are three 
linear independent equations for non-vertical incident 
15 waves and thus can be used to solve for the three 
unknowns i, a and y (and hence p) in a minimization 
scheme. Note that equations (61) and (62) should be 
normalized with slowness. 

No vertical derivatives are required so far, but if 
20 such information is available, similar equations can be 
derived for the corresponding combinations of vertical 
and temporal derivatives. In addition to this, 
expressions can also be obtained for an incoming S- 
wave, and these, two options enable us to improve the 
25 performance of the method, since the method is likely 
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to be more robust when larger parts of the data set can 
be used as input. In order to find more suitable events 
a measure for data-quality is required. A potential 
strategy is to use equation (63) and define an energy 
5 function as: 

- - [ - 

Once suitable events have been determined, the: 
inversion can be carried out over all these windows. 
10 The unknowns a and y should be equal for all intervals, 
whereas the angle of incidence is different for each 
event : 

7j«V.jr) = 2 s * j (64) • 

15 

where N is the total number of selected suitable 
windows and the E* are the energy functions using the 
k r / t interval. These are given by the squared differences 
of the left hand and right hand side expressions in 
20 equations (60-62) in case a time window is used with 
one incoming P-wave. 

Figure 7 is a flow chart showing the steps of 
estimating material properties according to preferred 
embodiments of the invention. In step 210, the seismic 
25 survey is designed using dense receiver groups as 

described and shown herein. In step 212, the seismic 
data are recorded. In step 214 the' spatial derivatives 
of the wavefield are estimated, preferably all centered 
at a single location using finite-difference and where 
30 desirable, the Lax-Wendroff corrections. In step 216 
the derivatives are used to obtain a system of 
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equations, or possibly a single equation. In step 218 
the results of the system equations are solved to find 
the most consistent material and/or wavefield 
properties and/or receiver coupling properties. 
5 Steps 220, 222, 224, and 226 show examples of 

possible uses of the results of step 218. In step 220, 
the property estimates are used to aid statics 
estimation. In step 222 the material properties are 
used to construct wavefield filters. In step 224 the 

10 receiver coupling properties are used to correct parts 
of, or the whole of the recorded wavefield. In step 
226 the wavefield properties are used to estimate 
subsurface reflection points . 

Figure 8 is a flow chart showing further 

15 detail of the preferred methods for estimating material 
properties, according to preferred embodiments of the 
invention. In step 300 the seismic survey is designed 
using dense receiver groups. The survey is designed as 
either a surface design, step 310, or as a volumetric 

2 0 design, step 33 0. In step 310, the entire receiver 

group is on the ground (or sea floor) surface. In step 
312, the seismic data are recorded. In step 314 the 
portion of the recorded wavefield that contains only a 
single in-coming arrival that is approximately planar 
25 is isolated. In step 316 the first order derivatives 
are estimated using the free surface condition to 
obtain dV/dz . In step 318 the estimates are inserted 
into plane-wave expressions to obtain a system of 
equations. In step 320 the resulting system of 

3 0 equations are solved to find the most consistent P and 

S velocities and arrival trajectory. 

In the case of volumetric seismic survey 
design, step 330, the receiver group includes at least 
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one receiver buried beneath the ground surface. In 
step 332 the data is recorded. Steps 342, 344, 346, 
and 3 48 depict the embodiment involving full wave 
equation inversion. In step 342, all first order, and 
5 all horizontal second and third order derivatives are 
estimated at a point P by finite difference ♦ In step 
344 the values of dV/dz and dV/dzz are re-centered at 
point P on the surface using- the Lax-Wendrof f 
correction. In step 346 all derivative estimates are 
10 inserted into the wave equation to obtain a system of 
equations, one for each point in time. In step 348, 
the resulting system of equations are solved to find 
the most consistent P and S. velocities. 

Steps 334, 33 6, 338, and 340 depict the 
15 embodiment involving vertical derivative inversion. In 
step 334 all. first order, and all horizontal second 
order derivatives are estimated by finite difference. 
In step 336, the values of dVz/dz are re-centered at 
point P on the surface using the Lax-Wendrof f 
20 correction. In step 33 8 the values are equated to 
expressions for dVz/dz given by the free surface 
condition, giving a system of equations. In step 340 
the resulting system of equations are solved to find 
the most consistent P and S velocities. 
25 While preferred embodiments of the invention 

have been described, the descriptions are merely 
illustrative and are not intended to limit the present 
invention. For example, while the preferred 
embodiments of the invention have been described 
3 0 primarily for use on the land surface, the invention is 
also applicable to receivers placed on and below the 
ocean floor. In the case of ocean bottom receivers, it 
is preferable to use stress conditions relevant for 
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fluid-solid boundaries rather than the free surface 
condition. Additipnally , the present invention is 
applicable to seismic measurements made in a borehole, 
known as borehole seismics. Although the examples 
5 described assume an essentially isotropic medium in the 
near surface region, the invention is also applicable 
to anisotropic media. In the case of anisotropic 
media, one may wish to increase the number geophones 
per group. 



10 
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CLAIMS 

What is claimed is: 

5 1. A method of estimating near- surface material 

properties in the vicinity of a locally dense group of 
seismic receivers comprising .the steps of: 

receiving data measured by the locally dense 
group of seismic receivers, the data representing 
10 earth motion caused by a seismic wavefield; 

estimating local derivatives of the wavefield 
such that the derivatives are centered at a single 
location in the vicinity of the receiver group; 
and 

15 using physical relationships between the 

estimated derivatives to estimate near- surface 
material properties in the vicinity of the 
receiver group. 

20 2. The method of claim 1 wherein the step of 

estimating the local derivatives includes centering the 
derivatives using a technique substantially similar to 
the Lax-Wendroff correction. ■ 

25 3. The method of. claim 1 wherein the physical 

relationships used to estimate material properties 
include the free surface condition. 



30 



4. The method of claim 1 wherein the physical 
relationships used to estimate material properties 
include wave equations. 
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5. The method of claim 1 wherein one at least 
one of the receivers in the receiver group is buried 
below the earth surface such that the receiver group 
encloses a volume. 

5 

6 . The method of claim 1 wherein the receivers 
in the locally dense group are spaced around 1 meter or 
less from each other. 

10 7. The method of claim 6 wherein the receivers 

in the group are spaced around 0.5 meters or less from 
each other. 

8 . The method of claim 7 wherein the receivers 
15 in the group are spaced around 0.25 meters or less form 
each other. 



9 . The method of claim 1 wherein the receivers 
in the group are spaced apart by a distance 
approximately equal to or less than one- fifth of the 
shortest wavelength of interest. 



20 



10. The method of claim 1 wherein the material 
properties estimated include those experienced by the 
25 incoming wavefield that are frequency and wave type 
dependant . 

11 - The method of claim 1 wherein' the receivers 
are located at or near the sea bottom, and wherein the 
30 physical relationships used to estimate material 

properties include fluid-solid boundary conditions. 
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12 . The method of claim 1 wherein the receivers 
are located in a borehole. 

13. The method of claim 1 wherein the physical 

5 relationships used to estimate material properties are 
derived from the physics of plane waves arriving at the 
receiver group. 

14. The method of claim 13 wherein the group of 
10 receivers does not include any buried receivers. 

15. An apparatus for estimating near-surface 
material properties in the vicinity of a locally dense 
group of seismic receivers comprising: 

15 a storage configured to store data measured 

by the locally dense group of seismic receivers, 
the stored data representing earth motion caused 
by a seismic wavefield; 

a derivative estimator adapted to process the 

2 0 stored data and estimate local derivatives of the 

wavefield such that the derivatives are centered 
at a single location in the vicinity of the 
receiver group; and 

a material property estimator adapted, to 

25 estimate physical relationships between the 

estimated derivatives to estimate near-surface 
material properties in the vicinity of the 
receiver group . 

30 16. The apparatus of claim 15 wherein the 

derivative estimator centers the derivatives using a 
technique substantially similar to the Lax-Wendrof f 
correction, and the physical relationships used to 
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estimate material properties include the free surface 
condition. 

17 . ] The apparatus of claim 16 wherein the 
5 receivers in the locally dense group are spaced around 
1 meter or less from each other. 

18 . The apparatus of claim 15 wherein the 
♦ receivers are located at or near, the sea bottom, and 
10 wherein the physical relationships used to estimate 
material properties include fluid- solid boundary 
conditions. 

19. The apparatus of claim 15 wherein the 
15 physical relationships used to estimate material 

properties are derived from the physics of plane waves 
arriving at the receiver group. 
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